A CTRW approach to normal and anomalous reaction— diffusion processes 
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We study the dynamics of a radioactive species flowing through a porous material, within the 
Continuous- Time Random Walk (CTRW) approach to the modelling of stochastic transport pro- 
cesses. Emphasis is given to the case where radioactive decay is coupled to anomalous diffusion in 
locally heterogeneous media, such as porous sediments or fractured rocks. In this framework, we 
derive the distribution of the number of jumps each particle can perform before a decay event. On 
the basis of the obtained results, we compute the moments of the cumulative particle distribution, 
which can be then used to quantify the overall displacement and spread of the contaminant species. 



I. INTRODUCTION 

The investigation of transport processes in inhomo- 
geneous geological formations has attracted intense re- 
search efforts, because of its relevance in the context of 
subsurface waste management and environmental reme- 
diation [E H, H, 0, [1] . In such complex physical systems, 
the spread of the transported quantity is often experi- 
mentally found to exhibit a non-linear growth with re- 
spect to time, of the kind (x 2 {i)) ~ f, 7 / 1. This 
scaling is actually the hallmark of the so-called anoma- 
lous diffusion, as opposed to Fickian (normal) diffusion, 
where 7 = 1 0, @ . 

The migration of contaminant particles through both 
homogeneous and heterogeneous materials has been suc- 
cessfully described within the Continuous-Time Random 
Walk (CTRW) scheme. In this stochastic model, the tra- 
jectory of each particle is represented as a series of ran- 
dom jumps separated by random waiting times, during 
which the walker stays at rest in the previously reached 
position 0, 0, IS & GH • F° r sake of simplicity, we adopt 
the common assumption that jumps and waiting times 
are independent of each other [3,[ll]]. The jump lengths 
are usually drawn from a Gaussian distribution with (fi- 
nite) variance a 2 , where a is a typical spatial scale de- 
pending on the traversed material, and mean [i 0, . 
A forward bias fi > is often used to model the con- 
tribution of an external advection field fill ]. I n the con- 
text of underground particle flow through porous sedi- 
ments or bedrock, the migrating plume is most frequently 
characterized by an anomalous spread (7 =/= 1), induced 
by the presence of, for instance, dead ends, stagnation 
and obstacles, which affect the particle dynamics at the 
microscopic scale [HI, [IH, [TH, EH [HI- These processes 
are mirrored in extremely long trapping times, which, 
within the CTRW formulation, are modelled by assum- 
ing that the waiting time distribution has a power-law 
decay w(t) ~ f- 1 - 01 with < a < 2 The 
broad distribution of spatial length scales which char- 
acterizes heterogeneous materials can result in a broad 
(power-law) distribution of characteristic time scales, so 
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that extreme events, i.e. anomalously long resting times, 
have a non negligible probability of being sampled. This 
phenomcnological picture is at the basis of the CTRW 
formulation 0, g M, [13] . 

In the case of independent jumps and waiting times, 
the general form of the CTRW transport equation for 
the normalized particle concentration P(x, t) can be ex- 
pressed as follows 

where the time convolution operator M, with a kernel 
M(t), takes into account possible non-Markovian (mem- 
ory) effects due to power-law waiting times (see Appendix 
lAl for details) [ll|. In particular, one-dimensional 
transport with a constant bias is subdiffusive (7 < 1) 
when < a < 1/2 and superdiffusive (7 > 1) when 
1/2 < a < 2, as shown, e.g., in [HGiJUO]. 

On the other hand, transport in a locally homogeneous 
material can be described by assuming that the asymp- 
totic decay of w(t) is sufficiently fast (as it is the case 
of an exponential distribution), so that the particles wait 
on average the same characteristic time between any suc- 
cessive jumps 0, H, 0]. In this case, Eq. (fT]) reduces 
to the well-known normal advection-diffusion equation 
0ll,[llj]. No te that the general formalism of CTRW can 
account also for a transition from anomalous to normal 
diffusion, by adopting for instance a truncated power-law 
distribution with an exponential cut-off for the waiting 
times: this behavior is often observed in contaminant 
migration (see, e.g., 0, HO, H [H ) . 

The theoretical framework of CTRW is well established 
and has been corroborated by a hug e amount of exper- 
imental evidences 0, [H Hzl El, |M H HI . However, 
due to the subtleties involved in the non-Markovian na- 
ture of the memory kernel 0, 12^ . |30| . much ingenuity 
has been necessary to couple reaction phenomena with 
anomalous diffusion [3J [H, [H, H3]- A comprehensive 
theoretical treatment, though, is still lacking: see, e.g., 
[3l| and references therein. 

In this paper, we consider the simple but significant 
case of a system composed of two diffusing species, say m 
and n, where m is unstable and decays through a nuclear 
reaction to n, which is stable. The decay is governed 
by a Poisson process with parameter A. This system 
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FIG. 1: The distribution V{N) (Eq. ([21]), solid line) is com- 
pared with Monte Carlo simulation (dots) for the following 
parameters: 10 simulated particles, to = 2 and r = 0.1. 




FIG. 2: The distribution V(N\T) (Eq. (J25J), solid line) is 
compared with Monte Carlo simulation (dots) for the follow- 
ing parameters: 10 5 simulated particles, T = 2 and r = 0.1. 



can characterize, e.g., the transport of a radioactive con- 
taminant species leaking from an underground repository 
and migrating through the surrounding geological forma- 
tions. In analogy with the well known normal reaction- 
advection-diffusion equations, it would be tempting to 
postulate a generalization of ([1]) with a decoupled struc- 
ture of the kind 

jfiPjfa t) = Mj]CjPj(x, t) t \P m (x, t), (2) 

where ICj = ijjd^/2 — (J,jd x is the transport operator and 
j = m,n . However, suitably extending the derivation 
of the CTRW scheme presented in [3l| it is possible to 
show that the concentrations of m and n obey to 

d 

OjPm(x,t) = M*JC m P m {x,t) - \P m (x,t) (3) 

and 

^-P n (x,t) = M n K n P n (x,t) + XP m (x,t), (4) 
at 

respectively, where the operator 

M* m = e- xt M m e xt (5) 

involves also reaction (A) terms [U IHj]. Thus, only 
the equation for the species n has a decoupled structure, 
where transport and reaction act independently. It can 
be shown that when Wj (t) is an exponential distribution 
the standard reaction-advection-diffusion equations are 
recovered, namely 

^Pjfat) = TjPjfat) T \P m (x,t), (6) 

where Tj = Djd% — Vjd x is the transport operator and 
Dj, vj are the diffusion coefficient and the velocity of 
each species j = m, n, respectively (3l| . 



We have implicitly assumed that particles m can still 
undergo a nuclear reaction when trapped in a stagnant 
region, and further that particles n once created have dif- 
ferent physical-chemical properties from m: these repre- 
sent reasonable hypotheses in the context of radionuclides 
migration. The concentration profiles corresponding to 
equations and ©, respectively, have been contrasted 
in [3l|: discrepancies are clearly visible, so that in prin- 
ciple it should be possible to select the proper model on 
the basis of available experimental data. Other possible 
implementations of reaction-diffusion phenomena within 
the CTRW formulation exist (see, e.g., [U [36j]), rely- 
ing on different physical assumptions and thus leading to 
different transport equations. 

Having this framework in mind, in the following we 
address the issue of computing the number of jumps a 
diffusing particle m can perform before decaying to n, 
and the corresponding overall displacement and spread 
of the radioactive species. In Section |TT] we outline the 
mathematical formalism. Then, in Sections IIIII and IIVI 
we discuss the case of normal and anomalous diffusion, 
respectively. Conclusions are finally drawn in Section [V] 



II. NUMBER OF JUMPS BEFORE DECAY 

Assume that the waiting times between consecutive 
jumps are sampled from independent and identically dis- 
tributed probability density functions (pdf's) w(t). Let 
w(u) — £{w(t)} denote the Laplace transform of w(t). 
Then, the distribution wjv(i) of t after jumps will be 
given by the TV- fold convolution of w(t) with itself: in 
the transformed space, we simply have wn(u) = w(u) N . 
Define 7 , (A|T) as the probability that a particle whose 
waiting times are distributed according to w(t) performs 
iV jumps before a final time T. The basic relation be- 
tween the counting process 7 5 (A f |r) and the pdf w(t) of 
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FIG. 3: The distribution V(N) (Eq. ([27]), solid line) is com- 
pared with Monte Carlo simulation (dots) for the following 
parameters: 10 J simulated particles, a = 0.5, tq — 4 and 
r = 10" 3 . 
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FIG. 4: The distribution V(N\T) (Eq. (|30]), solid line) is 
compared with Monte Carlo simulation (dots) for the follow- 
ing parameters: 10 s simulated particles, a — 0.5, T = 2 and 
r = 10" 3 . 



the waiting times between consecutive events is 

V(N\T) = W N (T)~W N+1 {T), (7) 

where Wn{T) is the cumulative distribution associ- 
ated to tujv(i), evaluated at T [37|. In Laplace space, 
V{N\u) = vr l (w N (u) - w N+1 (u)) = u^iwiu) 1 * - 
w(u) N+1 ). Therefore we have 

V(N\T) = CT X | i (w(u) N - w(u) N+1 ) | (T). (8) 

Let now f(T) = Ae _AT be the pdf of the radioactive 
decay events. Then, the probability that particles m per- 
form N jumps before decaying to n is 

f °° 

V{N) = I V(N\T)f(T)dT. (9) 
Jo 

Integrating once by parts we get 

/>oo 

V(N) = I e-^C- 1 {w(u) N } (T)dT + 

e-^C- 1 {w(u) N+1 }(T)dT. (10) 

Thus, interpreting each integral as a Laplace transform 
evaluated at u — A with respect to the internal argument 
C- 1 {w{u) N } (T), we finally have 



where 



V{N) = w{\) N -w(X) 



(11) 



Now, in order to characterize the displacement and the 
spread of the radioactive species m before decay, we are 
interested in computing the moments of the cumulative 
particle distribution P^{x), namely 



E[x r ] = \J x r P^(x)dx, 



(12) 



P m (x,t)dt 



(13) 



and the factor A is used to normalize the moments to the 
average radionuclide decay time. These quantities can 
be intuitively expressed in terms of the moments of the 
particles locations pdf after N jumps, pjv(x), averaged 
on the distribution V{N): 



E[x r ] = 22 ?W / x r PN {x)dx, (14) 

This can be understood as follows. First, note that, 
if Q m (x,t) satisfies Eq. ([T]) (without radioactive decay), 
then P m (x,t) — Q m (x,t)e~ xt satisfies Eq. © for the 
reactive species. Within the CTRW formalism, the con- 
centration Q m (x, t) can be expressed as 



JoW N (t') 



1 



f^' w(t")dt" ) dt' 



(15) 



where the quantity between square brackets corresponds 
to V(N\t) in Eq. © (see, e.g., [2l|). Then, integrating 
P m (x,t) over time (so to obtain the cumulative distri- 
bution P^(x)) and computing the r-th moment finally 
leads to expression (p~4]) . 

Assuming now that the single jump length has a Gaus- 
sian distribution with variance a 2 and mean /x, then 
Pn(x) is again a Gaussian distribution with variance Ncr 2 
and mean N/i. Therefore, the first and second moment 
of the cumulative particle distribution respectively read 



E[x 2 ] 



E[x] = fi (N) 
= o 2 (N)+n 2 (N 2 ), 



(16) 
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where brackets denote the average with respect to V{N). 
Finally, the radioactive species displacement is provided 
by the first moment E[x], whereas its spread can be 
expressed on the basis of the second centered moment 
S = E[x 2 ] - E[x} 2 [J . 

Furthermore, the link between P m (x,t) and Q m (x,t) 
allows the moments E[x r ] to be expressed as a function 
of the memory kernel M (f ) . Note indeed that Pf n (x) 
can be represented in terms of the Laplace transform of 
Q m (x,t), namely 



Pm{ x ) — Qm{ x ,X). 



(17) 



Then, it immediately follows that the moments of P^ (x) 
are given by the Laplace transforms of the moments of 
Qm{x, t). General expressions for multidimensional cases 
are provided, for instance, in [2ll |: in one dimension, we 
have 



E[x] = fj,X~ 1 M(X) 
E[x 2 ] = ( r 2 A- 1 M(A) + 2^ 2 A- 2 M(A) 2 . 

III. NORMAL DIFFUSION 



(18) 



We can now specialize this general formalism to the 
case of normal and anomalous diffusion. Within the 
CTRW approach, normal diffusion is usually modelled 
assuming that w(t) is an exponential distribution with 
mean r [H, El- In this case, the Laplace transform 
reads w(u) — (1 + ut) _1 , so that the kernel is simply 
M(u) = t _1 . Moreover, the convolution W]y(t) is known 
analytically and is given by the Gamma distribution (37| 



w N (t) 



N-l. 



t n t(n) 



whose Laplace transform reads 



wn(u) = w{u 



(1 + ut) 



N ' 



We can therefore obtain V(N): 



V(N) = w(X) N - w(X) 



\N+1 



t/tq 



(l + r/ro)^ 1 ' 



(19) 



(20) 



(21) 



where tq = A^ 1 . A numerical example is provided in 
Fig. [TJ where we compare Eq. (|2ip with Monte Carlo 
simulation. For each simulated particle, a random decay 
time T is first drawn from an exponential pdf with mean 
To. Then, the particle trajectory is followed until the 
cumulative waiting time (each contribution being drawn 
from an exponential pdf with mean r) exceeds T, and 
the number of performed jumps is recorded. Parameter 
values are provided in the figure caption. Finally, noting 
that EZokq k = 9/(1 - q? and ££o *V = q(l + 
q)/(l — q) , provided that \q\ < 1, we can compute the 
moments 



E[x] = ^ro 



E[x 2 



-TO 



(22) 



We assume that the time scale of transport is shorter 
than the time scale of decay (r -C To), hence the ap- 
proximation sign. In formula ([22]) . a 2 /2t is the diffusion 
coefficient D m and fi/r is the local particle velocity v m 
(induced by the forward bias /i) appearing in Eq. ([6]) 
hj. The same result could be obtained by resorting to 
expression (I18[) and substituting the specific functional 
form of M(u). 

When w(t) is an exponential pdf, the cumulative dis- 
tribution Wjv(t) is known exactly, namely 



W N (t) 



T(N) ' 



(23) 



where j(N,x) = J Q s N x e s ds is the (lower) incomplete 
Gamma function. Then, wc can explicitly compute 



V(N\T) = 



7 (iV,g) jjN + l^) 
T(N) T(N + 1) 



(24) 



This formula can be simplified by resorting to the prop- 
erties of the incomplete Gamma function, namely 7 (AT + 
l,q) = N"/(N : q) - q N e- q We thus get 



V{N\T) = 



N\ 



(25) 



which is a Poisson distribution with parameter T/r, 
as expected: V(N\T) is indeed a counting process for 
Markovian events whose average rate is t _1 , over a time 
interval T . A numerical example is provided in Fig. [21 
where we compare Eq. (|25[) with Monte Carlo simulation, 
which proceeds as in the previous case, provided that the 
random decay time is replaced by a fixed threshold T. 
Parameter values are given in the figure caption. 

These results can be extended to a broader class of 
distributions. It can be shown that any waiting time 
pdf with finite first moment would lead to an expan- 
sion of the kind w(u) ~ 1 — c\ut to the first order in 
u, i.e. sufficiently far from the source (ut < 1) @. The 
constant c\ > depends on the functional form of the 
pdf. To provide an example, for a Pareto distribution of 
the kind w(t) = ar a t~ l - a , with a > 1, we would have 
w(u) = 1 — c\ut + o(u a ), with ci = a /(a — 1). In order 
to recover normal diffusion, finiteness also of the second 
moment of the pdf w(t) is required in case of a non van- 
ishing bias /i, which therefore implies a > 2 (Lsl. IT9L [20| ■ 
Then, it follows that w(u) ~ 1/(1 + ciut) n and formu- 
las (|2~lj) and (|2"2")) . which have been derived for the expo- 
nential distribution, would remain asymptotically valid, 
provided that we replace r — > c\t. 



IV. ANOMALOUS DIFFUSION 

To illustrate the case of anomalous diffusion, a con- 
venient choice is assuming w(u) — 1/(1 + (itr) Q ), with 
< a < 1, so that w(t) ~ i _1_Q , for t — > 00, and the 
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kernel reads M(u) = u 1 a /r" [Hll]- The parameter r is 



a characteristic time constant. Then, 



wn(u) = w(u) 



N 



1 



(1 + {ur) a ) N 
and we can therefore easily compute V{N): 

(r/r r 



V(N) = w(X) N - w(\) 



\N+1 



(1 + (r/r )«) 



(26) 



(27) 



where tq = A -1 as before. A numerical example is pro- 
vided in Fig. [3J where we compare Eq. ([2T)l with Monte 
Carlo simulation for a — 0.5. The simulation proceeds 
similarly as in the case of normal diffusion, the waiting 
times being now drawn from a power-law pdf. Parameter 
values are provided in the figure caption. We finally get 
the moments 



E[x] = M 
E[x 2 ] ~ 2^r» + 2 



) 2 Tn 2 °. 



(28) 



Similarly as in the case of normal diffusion, we assume 
that transport occurs on a time scale shorter than the 
time scale of decay (r <C t ), hence the approximation 
sign. In formula (|28[). a 2 /2r a can be regarded as the 
generalized diffusion coefficient D* m — a 2 l /2r a and /i/r" 
as the generalized local particle velocity — [i m /T a 
implicitly appearing in Eq. Q 9|,lll|. This is true for the 
particular functional form of the memory kernel adopted 
here. The same result could be obtained by resorting to 
expression (|18[) and substituting the specific functional 
form of M{u). 

In this context, the long time behavior of the reactive 
species concentration P m {x, t) can be explicitly obtained. 
For the case of a vanishing bias (/i = 0), note that the 
contaminant concentration Q m (x, t) (without radioactive 
decay) can be expressed in closed form by means of the 
Fox's H function 



4-D* t° 



(I — a/2, a/2) 
(0,1) 



(29) 



provided that the solution is evaluated sufficiently far 
from the source [§, |3{|. The H function admits 
a computable representation as a series expansion, 
with an exponentially stretched decay \ogQ m (x,t) ~ 

- (\x\/t a l 2 ) 1,[1 ~ a/2) HIH. Then > the asymptotic prop- 
erties of P m (x,t) immediately follow from P m (x,t) = 

Qm(^i t)c 

In some specific cases, analytic results can be obtained 
for the distribution V(N\T). To provide an example, 
for the Levy-Smirnov pdf w(t) = (T/47r) 1 / 2 e _T/4t i _3/2 , 
which has a power-law decay with a = 0.5 [37| . the in- 
verse Laplace transform appearing in Eq. ([5]) can be 
explicitly evaluated, so that V(N\T) can be expressed in 
closed form as 

r (Nm ,J^jL).J-jL), ( 3„, 



where (p(x) — 2ir~ 1 / 2 e~ s ds is the error function. A 
numerical example is provided in Fig. 01 where we com- 
pare Eq. pop with Monte Carlo simulation. Parameter 
values are provided in the figure caption. In the general 
case, V(N\T) can be computed from definition @ with 
arbitrary accuracy by resorting to a numerical inverse 
Laplace transform algorithm [4(| ■ 

Similarly as for the case of normal diffusion, it can 
be shown that any pdf with power-law decay and infi- 
nite first moment would asymptotically lead to a Laplace 
transform of the kind w(u) = 1 — c q (ut) q + o(u), trun- 
cating the expansion to the first non constant term for 
iit « 1, i.e. evaluating the contaminant concentration 
sufficiently far from the source @ . The constant c Q > 
depends on the specific details of w(t): for the case of the 
Pareto pdf, for example, c a = T(l — a). The expression 
of w(u) can be regarded as the first order expansion of a 
pdf w(u) ~ 1/(1 + (wr) Q ). Therefore, formulas (f2"T|) and 
(f2"8")l would remain asymptotically valid, provided that we 
replace r a — * c a r a . 



V. CONCLUSIONS 

In this paper we have considered reaction-advection- 
diffusion processes within the CTRW framework, in both 
homogeneous and heterogeneous media, the latter giving 
rise to anomalous transport for the migrating species. We 
have analytically derived the distribution of the number 
of jumps that each particle can perform before under- 
going a reaction event. On the basis of this result, we 
have determined the moments of the cumulative particle 
concentration, which allow the overall displacement and 
spread of the reacting species to be quantified. Though 
we have focused on the case of radioactive contaminant 
particle transport, by virtue of its interest in the field of 
nuclear waste migration from underground repositories, 
the proposed framework could be applied to other phys- 
ical systems where the reaction term is linearly propor- 
tional to the concentration of the reacting species, such 
as first-order chemical reactions. 
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APPENDIX A: THE MEMORY KERNEL 

Let us briefly recall the definition of the Laplace trans- 
form: 

/>OC 

C{g{t)}{u)=g{u)= / e- ut g{t)dt. (Al) 
Jo 
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The convolution operator M. is defined as 

r+oo 

Mg(t) = / M(t - t')g(t')dt', (A2) 



where the kernel M(t) in the Laplace transformed space 
satisfies 



- , . w(u) 

M(u) = u 

y ' 1 - w(u) 



(A3) 



for a sufficiently well behaved function g(t) 0, [TTj] . It 
immediately follows that 



The properties of M. depend on the waiting times dis- 
tribution w(t). In the direct space, when w(t) has an 
algebraic decay, M(t) asymptotically behaves as a power- 
law kernel, accounting for long time correlations: these 
in turn induce non-Mar kovian (memory) effects. On the 
contrary, when w(t) is an exponential distribution the 
operator reduces to a constant, independent of time, so 
that the memory effects disappear, the transport pro- 
cess becomes Markovian and normal diffusion is recov- 
ered dEH]. 



£{Mg(t)} = M(u)g(u) 



(A4) 
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